{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Concept\n",
    "==============\n",
    "\n",
    "\n",
    "$$\n",
    "\\newcommand{\\sumN}{\\sum_{n = 1}^N} \n",
    "\\newcommand{\\sumn}{\\sum_n} \n",
    "\\newcommand{\\bx}{\\mathbf{x}} \n",
    "\\newcommand{\\bbeta}{\\boldsymbol{\\beta}} \n",
    "\\newcommand{\\btheta}{\\boldsymbol{\\theta}} \n",
    "\\newcommand{\\bbetahat}{\\boldsymbol{\\hat{\\beta}}} \n",
    "\\newcommand{\\bthetahat}{\\boldsymbol{\\hat{\\theta}}} \n",
    "\\newcommand{\\dadb}[2]{\\frac{\\partial #1}{\\partial #2}} \n",
    "\\newcommand{\\by}{\\mathbf{y}} \n",
    "\\newcommand{\\bX}{\\mathbf{X}}\n",
    "$$\n",
    "\n",
    "**<u>Model Structure</u>**\n",
    "\n",
    "*Linear regression* is a relatively simple method that is both extremely widely-used and a stepping stone for many more sophisticated methods. This makes it a natural algorithm to study first.\n",
    "\n",
    "In linear regression, the outcome variable $y$ is assumed to follow a linear function of one or more input variables, $x_1, \\dots, x_p$, plus some error. Specifically, we assume the model for the $n^\\text{th}$ observation in our sample is of the form \n",
    "\n",
    "\n",
    "$$\n",
    "y_n = \\beta_0 + \\beta_1 x_{n1} + \\dots + \\beta_px_{np} + \\epsilon_n. \n",
    "$$\n",
    "\n",
    "\n",
    "Here $\\beta_0$ is the intercept term, $\\beta_1$ through $\\beta_p$ are the coefficients on our input variables, and $\\epsilon$ is an error term that represents the difference between the true $y$ value and the linear function of the predictors. Note that the terms with an $n$ in the subscript differ between observations while the terms without, namely the $\\beta\\text{s}$, do not.\n",
    "\n",
    "It is often easier to write our predictors as a vector. Let $\\bx_n = \\begin{bmatrix} 1 & x_{n1} & \\dots & x_{np} \\end{bmatrix}^\\top$ (we add the leading 1 to correspond to $\\beta_0$). Then we can equivalently write, \n",
    "\n",
    "\n",
    "$$\n",
    "y_n = \\bbeta^\\top \\bx_n + \\epsilon_n.\n",
    "$$\n",
    "\n",
    "\n",
    "Below is an example of a linear regression model with a single input variable. The input variable is generated randomly and the outcome variable is generated as a linear combination of that input variable plus an error term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEcCAYAAAAhoQi5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU1fnH8c8DYS8qoqAFKa6IC+C+IiKoaOv2U6tYURSqoohWERVXpFpQcKug4l5EqVoVl4pAXQBXlirKVoIiW2UrIBASEnJ+f5w7YRiyTJKZuTOZ7/v1yiuZmXvnPrmZ3Oeec899jjnnEBGR7FQr7ABERCQ8SgIiIllMSUBEJIspCYiIZDElARGRLKYkICKSxZQEqsDMepqZM7P9ylnmRTNblMKwEsLMWge/W+Rrq5n9bGZjzGyvsONLBTNbZGYvhrDd48zsKzPbFOz7DqmOId2Y2fFmVhzsj5wUbvcAM3vMzGaZ2UYz+6+ZvWNm7VMVQ6qkbKdmocHAY2EHUQ1/Ad4B6gLHAvcAbc3sGOdcYaiRJd95wC8hbPc5YDNwFpAH/CeEGNKGmdUBngZWAHukePOnAZ2Bl4CZwC7AAOArMzvBOTcjxfEkjZJAkjjnFoYdQ1mCf64iV/6dgj84574Mfp4crPNn4Ajgy7JXSywzqw2Yc64oVdt0zv07VduKMLNaQBvgfufcRwl4PwPqOOe2VDu48NwCGPA8MLAqb2BmrYEfgc7OuU8qsepYYET0/4iZfQQsAm4ALqtKPOlI3UFJEtsdFNXNcrWZ3Rc0L9eZ2btm1rKU9f9oZt+aWb6ZrTaz58xs15hl+prZF2b2v+C9vjSz38YsE9nutWb2oJktBwrwZzaVMTP43irm/Rua2VAz+9HMtgTf7wgOatHLHW5mU4LfZ4mZDTSzQWbmYpZzZna/md1mZj8CW4BDg9d2M7MnzWyZmRWY2Twzuypm/T3M7CUzWx4s818ze8/MmgWv55jZYDNbGLVvp5rZiVHvsUN3kJkdbWaTgq6BTWb2LzM7OmaZF81sqZkdFvyueWa2wMyuKW/HmllPYCv+//GuYB8sinr90pjPwmgz2zPmPRaZ2ctmdqWZzQv223afhahlvzOzt0p5/uRg26eXF28qmNm+wB3AtUDKW57OudWxJ0nOufX41lmLyHNm1jvYZ+dGPVfbzCYHn7HGqYu6ipxz+qrkF9ATcMB+5SzzIrAo6nHrYJ1FwCvAGcDlwGrg05h1h+A/+MPxzdIrgGXAV0DtqOWGAb2ALsDpwBPBNs4oZbvLgLeB3wHnAA3KiDuyfO+Y5/sEzx8R9VwOMAVYA9wYxHEHkA8Mj1puN2AtMBv4PXAu8Cnwk/8IbredSKxTgPOBbkBzYCdgPrAY+CPQFXgIf/C8Pmr9ifh/1D8AJwEXAk8BrYPX7wA24s/mOuG7XgYBZ0e9xyLgxajH7fDdNDOAC4K4pgXPtY/5m/8CzAWuBk4N/tYOfyZa1mdld+CEYLln8d1vhwWvXRU8PxY4E+gNrAx+x1/FxLwM+B7oHvwt9i1je5ED669jnn8V+AHf8irv858Tx1ft8t4jjv+xCcBLwc/3Bvsgpwrv0zpY9+QE/N/vCmwCnoh5/jX8/3GLqHgLgWOqu81UfIUeQCZ+Ub0kEHvA7x88/+uo5bYCd8csFzlInFvG9moF/3wTgHGlbHdmRf/cMctfFbxfQ+AUYCnwRsyyPYJlT4p5/g78mWiz4PEDweOWUcs0wPf1uph1HbCcmCQF3IVPLvvHPP9M8A+YEzzeCPQr5/d7D3izgn2wiO2TwBvAOmCXqOd2Av4X/V7B33y7Az5QL4hvVAXbzAnWvTfqudrBPvo4ZtkTg2X7xcScB+wRx9+4MT5Z3RX13G74FuJtcX4+KvpaVFEc5Wzj0mDfRj4/9xJnEoj6P4h87Rus2yXm+VpViGtMsI/3i3l+F/wJzcf4E4si4Paq/v6p/lJ3UOq9H/P4u+B7pJvlVPwHeUzQdZFjflTEV/h/3JMiK5rZEUFXxwr8B68wWL9NKdt92wWf2Dg9HbzfJuBf+IPRpTHLdMN/+D+PiXUCUAd/Rkvw/Qvn3NLIis65zey4LyLGB6/Hbusr4MeYbX0INAUOCpabBtxiZjeY2aFB33i0acCZQZfTiWZWN459cRLwnnNuXVT8v+AvnHeKWTbPOfdx1HIFwAJiutHi1AZohj/4lHDOTcXv99htf+mc+7miN3XObQBeBnpHddtdge9/f6GC1ZcDR8XxdVZ5b2JmtaL/juav/RB0eQ4HBjrnVlb0u5TibvznNvKVGzw/Keb5uyvzpmZ2O3AJ0Nc5lxv9WvC5uAToiP88TgGGViH2UOjCcOr9L+ZxQfC9fvC9WfA9l9I1BTA/XPNfwBzgenw3SRF+VFLbUtb7byXj/DMwDt8SuBDoC4wEroxaphnwG8rus20afN8T300Ra0UZ65UWazNgvzi2dRF+JNMA4FHgv2b2FPBn51wxvlWSj09oA4GNZvYGcItzbnUZ771rGTH9DDSJeW5tKcsVsO3vWxmRa0BlbXvXmOcq8zceie/iO9PM3se3/N5yzpX1NwHAObfFzL6J4/0rOuG4G/93ivgUOBn/uVsBvGZmketWkX23s5nlO+c2lfO+o/CtvYg98cn6Gnx3XsTyCuIrEVzTeQC40zn3fBmLfYnvrjwIeCz4rGUEJYH0syb4fhqlH1Air3cDdgZ+H32GbWYNy3jfytYM/8k5Nz34eXJwgesKM3vKOfd1VCw/4vv5S7Mo+P5ftiW3aM0rEesafF/4DWWsMx8gOHu8DrjOzNrgr7sMAlYBTzo/vHUoMNTM9sBfI3kYn+wuKuO9/0fpQxT3YMeknkiR9y5r29Njnov7b+yc+97MpuCvXeTjE+zVFa1n20bbVOQnfNdRWWIP1huC7wfhBwKs2WEN3602Dn9NqVTOueVEHeCDeAHmR32e42ZmPfAJc7hz7v5yFr0H2B+YBTxiZh87fyE57SkJpJ+JQDHQyjk3sZzlIgf7kjNjMzsAf+1gaalrVM9t+IP9PWwbdTIef5F0o3NuXjnrfgn0N7OWkYRlZg0oY/RKGcYTtHji7SZwzs0HBgZncoeU8vrPwLNmdmZpr0f5FPitmTUOulIIkuJZwCeV+B0qaz7+rPhi/D0EBNs+Ht8CG17N9x+J7xZqAvzHxTc0NdIdVJGC8l6MPVhHuZEdR671xCfzrpTdekw4MzsP3z32rHOufznLdcS3Km8H/g58CzyJ7yJKe0oC1dPNzGL7YNdXcPAul3NuoZkNBZ4IzmQ/xZ+p7YXv73826HOehO/++ZuZDcc3ewfhu4USfq3HOfezmY3AH8yPcP5mmTH4vuR/BTF8i7+5bF/gbPxF7Dz8mXYf4EMzG4Q/QNwUfI/37PUR/Jn6FDN7BH+AbAQcCHR0zp1jZjvj98sYYB4+QZ6DP8hNADCzcUGcM/EtrcPwraqny9n2YHyL4V/B38YBt+IT8X1xxl9pzrmtZnY38LSZvYw/YLcA7sdfZ6io/74i/8B3mZ0A3BxnTFvYsQWSMM65HbqazOzk4MdPXYruFzGzk/CjpWYBL5rZsVEvF7jgXhIza4L/vH0MDHPOOfPDll8zsw+dcy+lIt7qUBKonr+W8txsyj+rrJBzbqCZzSXo1sAfdJbgrwEsCJaZbWZ/wB+E3gEW4s/Wu+H7VpNhCL7L4G7gHOdcYTCm/DZ8n/Le+AvJC/EXfbcEsa42sy7A48Df8E39p/AjUuK66cY5tz44A74bfwBugR+xMx9/MAOfLGfih5D+Bt+img/8wTk3LlhmMv4ax3X4g/hi4EH8gbWsbc8KDkT34+8gNXzrppNz7tt44q8q59woM8vD3zg1Dj/66Z/AAOfcxmq+d2GQFC/H/16yzSn4kV2HAZ/FvBbd1TUKP9LtssjAC+fc62b2HP5E7rPYC8npxio3YEQkMYLRIDOB1c65LmHHk42C0VW5wBTnXI+w45FwqCUgKWFmg/EHnJ/wI3l642/COjPMuLKRme2Eb61egu9mrO61BclgSgKSKg7flfPr4OdZ+GsGH4QaVXY6HN+HvRK4obR+eMke6g4SEcliumNYRCSLZVR3ULdu3dz48ePDDkNEJNPElk8pkVEtgdWry7qrX0REqiKjkoCIiCSWkoCISBZTEhARyWJKAiIiWUxJQEQkiykJiIhksYy6T0BEJFvkFRSxoaCItXlbaNKwLo3r5dCwXuIP2UoCIiJpJq+giElzV9D/9Vls2VpM3dq1GHZhO7q2bZ7wRKDuIBGRNLOhoKgkAQBs2VrMLW/MYkNB4ufUURIQEUkza/O2lCSAiIKiYtblFZaxRtUpCYiIpJkmDetSt/b2h+d6ObXYpWGdhG9LSUBEJM00rpfDsAvbUS/HH6Lr5dTioQva0zgJF4Yzaj6BI4880k2fnrQ5rkVE0kZkdNC6vEJ2aVinuqODyqwiqtFBIiJpqGFw0G++U/2kbkfdQSIiWUxJQEQkiykJiIhkMSUBEZEspgvDIiJpKhX1g5QERETSUKrqB6k7SEQkDUXqBzVbs5yDVyxMWv0gJQERkTS0bu0GrpryCpOeu5b7PxwBziWlfpCSgIhIupk0iX27nkD/KS8zad+juebcgWCWlPpBuiYgIpIuli+Hm26Cv/+dWvvuxxcjXubm5btSUFSctPpBSgIiImErKoK//hXuuQe2bIFBg6g1YADtLYfJiasfVColARGRMH32GVx7LcyaBWec4ZPBvvsC0BCSXj9I1wRERMKwahX06gUnnghr18Kbb8L775ckgFQJNQmYWTczm29muWZ2W5ixiIikRHExjBoFbdrA3/4Gt94Kc+fCeeeBlVnxOWlCSwJmVhsYAZwBHAR0N7ODwopHRCTpZs6E44+Hq6+Gdu3g229hyBBo1Ci0kMJsCRwN5DrnfnDObQHGAueEGI+ISHKsXw/9+sFRR8GPP8Lo0fDxx3BQ+Oe9YSaBFsCSqMdLg+e2Y2ZXmdl0M5u+atWqlAUnIlJtzsGYMb7rZ8QI6NMH5s+HSy8NpeunNGEmgdL2wA5zXTrnRjnnjnTOHbn77runICwRkQSYOxe6dPEH/Fat4Ouv4YknYJddwo5sO2EmgaXAXlGPWwLLQ4pFRCQxNm2C227zff7ffANPPQVffAFHHBF2ZKUK8z6BacD+ZrY3sAy4GLgkxHhERKrOORg3Dm64ARYvhp49YehQaNYs7MjKFVoScM4VmVlf4EOgNvC8c252WPGIiFTZDz/4C7/vvw+HHAKTJ0PHjmFHFZdQ7xh2zv0T+GeYMYiIVFlBATz4IDzwAOTkwLBhPhnUSWyRt2RS2QgRkaqYOBGuuw4WLIALLoBHHoGWLcOOqtJUNkJEkiqvoIgVv+Qz7+dfWPFLPnkJnhQl5ZYtg4sugtNO89cBxo+H11/PyAQAagmISBKlaorElCgs3Fbps7AQBg2CAQOgfvKKu6WCWgIikjSRKRK3bC0GSNoUiUn32Wd+iOfNN8NJJ8Hs2XD33RmfAEBJQESSaG3elpIEEJGMKRKTZtUquPJKX+lz3Tpf6fO991Je6TOZlAREJGmaNKxL3drbH2aSMUViwhUXw9NP+3IPo0f7bp8QK30mk5KAiCRN43o5DLuwHfVy/KEmWVMkJtTMmXDccXDNNdC+va/0OXRoqJU+kymN/xIikuka1suha9vmTB7QOalTJCbEunVw110wciTsvrtvAfzhDzXuzD9WGv4lRKQmaRgc9JM5RWK1OAevvOIv+q5a5ad6HDw47Qq9JYuSgIhkr7lz/UH/k0/g6KN92Yc0LfSWLLomICLZZ9MmuP32bbN7pXmlz2RSS0BEskdspc8rrvAXfbN4rhIlARHJDrGVPqdM8eP/s5y6g0SkZiso8Bd6Dz4YPv0Uhg/3w0CVAAC1BESkJouu9HnhhfDwwxlb6C1Z1BIQkZqntEqfr72mBFAKJQERqTkKC/3Z/oEHwjvv+Eqf330Hp58edmRpS91BIlIzTJ3qx/x/9x2ceaYv+7zPPpV6i7yCIjYUFLE2bwtNGtZN37ubE6hm/3YiUvOtWgW33govvAB77QVvvQXnnFPpcg81au6DSlB3kIhkprIqfZ57bpXq/dSYuQ8qSUlARDLPjBkJr/SZ8XMfVJGSgIhkjnXr4PrrfZ2fn36Cl1+Gjz6Cgw6q9ltn7NwH1aQkICLpzzl/wD/wQF/q+dprYd68hJZ6zsi5DxKgZv92IpL55sxh6zV9qD1lMpsPO4K8196iwTFHJfxibUbNfZBAagmISHratAluuw3Xvj1bv/2Wu864noNOvYfjPlzHpLkryEvCBdvIvAdt9mhM853q1/gEAGoJiEi6iVT67NcPliwh/9LLOLnpGayo3xjYNmpn8oDOWXGQTja1BEQkffzwA5x1lp/QfZddYOpUfnroryUJICIbRu2kipKAiISvtEqfM2bACSdk7aidVFESEJFwTZgAhx4Kd98NZ5/tR/3cdBPU8Qf5bB21kyraiyISjmXL/MH+tddg//3hww991c8Y2TpqJ1VC2Ytm9hBwFrAFWAhc4ZxbF0YsIpJihYW+uNs990BRke8GuuUWqFevzFUaBgf95jvVT2Gg2SGs7qCJwCHOuXbAf4DbQ4pDRFJp6lQ/mfvNN0OnTjB7Ntx5Z7kJQJIrlCTgnJvgnIsM8v0S0EwPIjXZqlV+UveOHWH9enj7bXj33UqXepbES4cLw1cCH4QdhIgkwdat2yp9vvwy3HYbzJlTpVLPkhxJuyZgZpOAPUp56Q7n3LhgmTuAImBMOe9zFXAVQKtWrZIQqYgkxYwZvsbP11/DySfDiBEJKfQmiWXOuXA2bHY5cA3QxTmXF886Rx55pJs+fXpyAxOR6lm3zvfzP/kk7L67n+6xe3ed+YerzJ0f1uigbsCtQKd4E4CIpDnnYMwYf9F39WrfChg82N/5K2krrIG2TwD1gInmzw6+dM5dE1IsIlJdc+bAddfBJ5/4Wv8ffACHHx52VBKHUJKAc26/MLYrIgm2aZM/2x8+HBo39heBe/eGWukw5kTioVvuRKTyYip9cuWVMGSIvwYgGUXpWkQqp5RKnzz3nBJAhlISEJH4xFb6fPhhmDkTTjgh7MikGtQdJCIVmzAB+vaFBQvgoov8NYAWLcKOShJALQERKduyZfD738Ppp/vHH34IY8cqAdQgSgIisqPCQt/dc+CBvsbPfffBd9+VWupZMpu6g0Rke1On+hu9vvsOzjzTl31WobcaSy0BEfFKq/T53nsZnwDyCopY8Us+837+hRW/5JNXUFTxSllELQGRbLd1KzzzDAwcCBs2+Eqfd94JjRqFHVm15RUUMWnuCvq/PostW4upW7sWwy5sR9e2zTUzWUAtAZFsNmMGHHcc9OkD7dvDt9/CX/6S0AQQ5pn4hoKikgQAsGVrMbe8MYsNag2UUCoUyXB5BUVsKChibd4WmjSsG9/8u5FKnyNHQrNmMHo0/OEPCa/0GfaZ+Nq8LSUJIKKgqJh1eYWaqjKgloBIBoscZDsO/Zhuj06h49CPmTR3Rdln2875yV3atPGlnvv2hXnz4NJLk1LqOewz8SYN61K39vaHuXo5tdilYZ2UbD8TKAmIZLBKHWRnz4bOnaFHD2jdGqZNg8cfT2qp5/LOxFOhcb0chl3Yjno5/lBXL6cWD13Qnsa6HlBCe0Ikg8XV3bFxoy/38PDDvtLnqFHQq1dKKn1GzsSjY0zlmXjDejl0bducyQM6sy6vkF0a1omvuyyLqCUgksHK7e5wDt56y0/p+OCDcNllMH8+/PGPKSv1nA5n4g3r5dB8p/q02aMxzXeqrwQQI7TpJatC00uKbC9yTeCWN2ZRUFRccpA9te4GGvT/k5/c5dBDff9/SIXeIheudSYeqjIv+CgJiGS47Q6ytbbS5IlHqfvgEKhTx5d7uP56yNFBN8ul1xzDIpI4DYMz6+ZfTvZTPObm+qJvDz+sQm9SIV0TEMl0S5duq/Rp5ss+//3vSgASFyUBkUwVqfTZtq2v9Dl4sC/6duqpYUcmGUTdQSKZKLrS529/6yt97r132FFJBlJLQCSTlFbp8913lQCkyipMAmZ2mZmtMbN6Mc+PMbN3kheaiJTYuhWeftqXexgzxlf6nDMHzjknKeUeJHvE0xJ4PVjunMgTZrYzcB7wXJLiEpGIGTPg+OPhmmuSVulTsleFScA5txkYA1wZ9fQlwC/A+0mKS0TWrfMF3o4+Gn76ybcAPvrIXwgWSZB4rwk8A5xqZi2Dx1cCLznnVJRbJNGc86WdI5U+r73WV/q85BJ1/UjCxTU6yDn3rZnNBHqa2dvAkcClSY1MJBvNnu0P+pMnwzHH+LIPhx8edlRSg1VmdNAzQE+gN/CZc25+UiISyUYbN8Ktt0KHDvD9977S5+efKwFI0lUmCbwK7AH0QReERRLDOXjzTd/P/+CDcPnlKa/0Kdkt7k+Zc24D8BqwJfguItWxcCH87ndw/vmw667+BrBnn4Xddgs7MskilT3V2BMY65zblIxgRLJCfr6v7nnwwb7v/5FH/DDQkEo9S3aL68Kwme0KdAVOA9onMgAz6w88BOzunFudyPcWiValCdkT7cMP/bDP3Fy46CIYPlyF3iRU8f4HzAR2BQY6575P1MbNbC/gVGBxot5TpDSRyVci8/HWrV2LYRe2o2vb5qlJBEuXwp/+BG+8Afvv7yt9qtCbpIG4uoOcc62dczs554YmePuPAAOAzJnZRjJSpSZkT6TCQn+237YtvPce/PnPqvQpaSW0KqJmdjawLLgHobzlrgKuAmjVqlWKopOaJq4J2RNt6lTo08cP+VSlT0lTSU0CZjYJP6w01h3AQPw1hnI550YBo8BPL5nQACVrRCZkj04EJROyJ9rKlX7M/4svQqtWvtLn2Wfrbl9JS0kdiOyc6+qcOyT2C/gB2Bv41swWAS2BmWZWWsIQqbbG9XIYdmE76uX4j3xkQvbGibwesHUrPPWUKn1KRkmLieaDRHBkRaODNNG8VMd2E7I3rJPY0UEzZviun2nToHNnGDFChd4knWiieZGSCdkTeQ1g3Tq4804YORKaNYOXX1ahN8koaZEEnHOtw45BpFKc8wf8/v1h9Wo/9n/wYNh556RsLi3ucZAaSZ8ikcpKcaXP0O9xkBpNFapE4hVSpc/Q7nGQrKAkIFKRkCt9lnePg0h1KQmIlCc319/oFan0+dlnKa/0GbnHIVrS7nGQrKMkIFKa/HwYNAgOOQSmTNlW6fP441MeSkrucZCslRb3CcRL9wlISowf70f7LFyYNpU+k3qPg2QD3ScgUqHoSp8HHAATJ0LXrmFHBSTpHgcR1B0k4it9DhsGBx7oK30OHgyzZqVNAhBJJrUEJLtNmeLLPcye7ad6fPxxVfqUrKKWgGSnlSuhZ0846SQ//n/cOHj3XSUAyTpKApJdoit9vvIK3H67bwWcfXbYkYmEQt1Bkj1U6VNkB2oJSM23bh1cdx3uqKPY+tNilj/xDCveep+8ffYPOzKR0CkJSM3lHIweDW3a4J56ih8v6slRlz7B8Uv2pOODnzBp7gryVH9HspySgNRMs2fDySfDZZfB3nvzv08+o9s+F/K/Og0AFWETiVASkJpl40YYMGCHSp+r9j9IRdhESqELw1IzRCp93nijv/O3Vy8YMqSk0FtKJ5oXySBqCUjmW7gQzjwTLrgAmjYttdKnirCJlE4F5CRz5efD0KHwl79A3bq+3MN110FO6Qd2FWGTLKYCclLDRFf6vPhiX+nz178udxUVYRPZkbqDpNryCopY8Us+837+hRW/5Cd32OWSJb7b54wzoHZtX+nz1VcrTAAiUjq1BKRaUjYJemEhPPYY3HuvL/1w//1w881Qr17itiGShdQSkGpJySToU6bAYYfBLbf4cg9z5sDAgUoAIgmgJCDVktRJ0GMqfRa88SYrxrzOvAZNk9/tJJIl1B0k1ZKU8fdbt/qbvAYOhE2b4Pbbyet/K5MWb6T/0I+T2+0kkmXUEpBqSfj4+xkz4Ljj4Npr4fDD/QxfDzzAhpx6ye92EslCOoWSamlYL4eubZszeUDn6o2/X7sW7rwTnnwSmjeHMWOge3cwP7y5vG4nDfkUqTolAam2ao2/dw5efhn694fVq+H66+G++2DnnbdbTGUfRJJD3UE1UErH7VdHTKVPpk/3w0BjEgCo7INIsqhsRA2TsnH71bFxoz/bf+QR2GknX/rhyiuhVvnnJCr7IFJlZZaNUEughknJuP2qcg7+8Q8/peNDD8Hll8P8+dC7d4UJACjpcmqzR2Oa71RfCUAkAUJLAmZ2vZnNN7PZZvZgWHHUNEkdt18dubkVVvqsjIzp8hJJc6GcSplZZ+AcoJ1zrsDMmoURR02UdhdQYyt9PvaYH/5ZRqXPeGREl5dIhgirJdAHGOKcKwBwzq0MKY4aJ60uoI4fD4cc4uv9nHcezJsH/fpVKwFAmnd5iWSYsE6bDgA6mtn9QD7Q3zk3rbQFzewq4CqAVq1apS7CDJWwcfvVsWQJ/OlPvv+/TRuYNAm6dEnY2+ueAZHESdqRwcwmAXuU8tIdwXabAMcCRwGvmdk+rpShSs65UcAo8KODkhVvTRJa3fzCQnj0URg0CIqL4c9/9uP/E1zoLe26vEQyWNKSgHOua1mvmVkf4M3goP+1mRUDuwGrkhWPJNnkyb6vf/ZsOOssePxxaN06KZuKdHnd8sYsCoqKdc+ASDWE9V/zNnAK8ImZHQDUBVaHFItUx8qVvsTz3/4Gv/kNjBsHZ5+d1E1Wp8srcq/B2rwtNGlYV/caSNYL69P/PPC8mX0PbAEuL60rSNLY1q3w9NNwxx2+0ufAgf7nhg1TsvmqdHlpVJHIjkL55DvntgCXhrFtSYBp06BPH1/x85RTYMQIOPDAUEOK5wy/rFFFkwd0VhKQrKVPvsRv7Vp/tv/UU77S56uvwkUXlVT6DEu8Z/gaVSSyI5WNkIo55/v827TxXUD9+vkx/xdfHHoCgPjvG4iMKoqmUUWS7ZQEpHzffw+dOvk6P/vu67uAHn201EqfYYm3VEZa3Ugnkib06ZfSxVb6fOaZuCp9VkaiRurEe99AWtxIJ5Jm9OmX7TkHbySwhtIAAA4sSURBVL4JN94IS5dCr14wZEiVC72VJZEjdSpz30BoN9KJpCnNJyDb5Ob6mb3Gj4f27WHkSDj++KRsasUv+XQMJo2PqJdTi8kDOlfpAK25BkTKVebFO/2XyI6VPh99FK67rtqF3sqT6JE6OsMXqRolgWw3fjz07QsLF/qJ3YcNg1//OumbVf0fkfSg0UHZaskSP8HLGWf4M/5Jk+CVV1KSAEAjdUTSha4JZJvYSp933gk335zwSp/xUD++SMromoCQ0kqf8VA/vkj41B2UDVas8Dd7derkx/+PGwfvvBNqAqhJNN+xZDK1BGqyrVth1Chf4XPTJrj9dt/9k6JKn9lAlUkl06klUFNNmwbHHuu7f444AmbNggceUAJIMM13LJlOSaCmWbvWH/iPOcbf8fvKKzBxYuilnmuqeOsWiaQrJYGaoqxKn927p0Wlz5pKlUkl0ykJ1ASzZ6d9pc+aSvc7SKbTfQKZbONGP97/0Ud9pc8HH4QrrkhopU+pmO53kAyg+wRqlNhKn717+7o/Ca70KfHR/Q6SyXTKmGlyc+HMM33Jh6ZN4fPPfa1/JQARqQIlgUyRnw/33guHHAKffQaPPQbTp8Nxx4UdmYhkMHUHZYIPPvB1/hcu9PP6Dh+eskJvIlKzqSWQzpYsgfPP990/kUqfr76qBCAiCaOWQDqKrfR5//2hVfoUqa7CwkKWLl1Kfn5+2KHUePXr16dly5bUqRP/fSpKAukmutLn2Wf7vn8VepMMtnTpUho3bkzr1q0x3biYNM451qxZw9KlS9l7773jXk/dQelixQq47DJ/09emTb7K57hxSgCS8fLz82natKkSQJKZGU2bNq10i0tJIGxbt/oJ3du0gbFj4Y47ttX7F6khlABSoyr7Wd1BYZo2Dfr08WUeunSBJ55QoTcRSSm1BMIQXelz+XI/4keVPkWSYs2aNXTo0IEOHTqwxx570KJFi5LHW7ZsSdh2Jk2ahJnx0ksvlTw3bdo0zIxHH3007vfJzc2lQ4cO1V4mXmoJpFKk0uctt8CaNb7S5333+bo/IpIUTZs25ZtvvgHg3nvv5Ve/+hX9+/ffbhnnHM45alWz7tahhx7K2LFjufzyywEYO3Ys7du3r9Z7JltWJIFIga+1eVto0rBuOAW+vv/en/1PmeLv8p0wARKUyUUyxo03QnBATpgOHfyQ6krKzc3l3HPP5cQTT+Srr77i7bffpn379qxbtw7wB/BJkybx7LPPsmLFCvr06cPixYupVasWjz/+OMcee+wO77nPPvuwatUqVq9eTdOmTZk4cSJnnHFGyeszZ86kT58+bN68mf3335/nn3+enXfemWnTptGrVy8aNWrECSecULJ8UVERAwYMYOrUqeTn59OvXz969+5dhZ1UttC6g8ysg5l9aWbfmNl0Mzs6GduJTP/XcejHdHt0Ch2HfsykuStSNw/sxo3+zL9DB3/B99lnYepUJQCRNDBnzhx69erFv//9b1q0aFHmcv369WPAgAFMnz6d1157rdwD8fnnn88bb7zB5MmTOeaYY7Ybs3/ppZcyfPhwZs2aRZs2bRg8eDAAPXv25Mknn+SLL75g69atJcuPGjWKZs2a8fXXXzNt2jRGjBjB4sWLE/CbbxNmS+BBYJBz7gMzOzN4fHKiN1LW9H+TB3RObmsgUunzhhtg2TJV+hSBKp2xJ9O+++7LUUcdVeFykyZNYv78+SWP165dy+bNm2nQoMEOy1500UX06NGDAw44gO7du/PRRx8B/tpEfn4+J554IgCXX345PXr0YPXq1WzevLmkBdCjRw8+/vhjACZMmMDcuXMZO3YsAOvXr2fBggX85je/qd4vHiXMJOCASGf4zsDyZGykvOn/klb6NzfX1/oZPx7at4fXX1ehN5E01KhRo5Kfa9WqRfT8KtHj7Z1zfP3119StW7fC92zRogXOOT799FNGjhxZkgTKm7ulrKGdzjlGjhxJly5dtns+Nze3wjjiFebooBuBh8xsCTAMuL20hczsqqC7aPqqVasqvZGUTv+3ebMqfYpkqFq1atGkSRMWLFhAcXExb731VslrXbt2ZcSIESWPv6ngusbgwYMZOnTodhead9ttNxo0aMDnn38OwOjRo+nUqRO77bYb9evX54svvgBgzJgxJeucfvrpjBw5kqIi3309f/58Nm/eXP1fNkpSWwJmNgnYo5SX7gC6AH9yzv3DzH4PPAd0jV3QOTcKGAV+ZrHKxhCZ/u+WN2ZRUFScvOn/oit9du/uK33uuWdityEiSTV06FC6detGq1atOOiggygoKABgxIgR9OnThxdeeIGioiI6d+68XVKIFenyiTV69OiSC8P77bcfL7zwAgAvvPACvXv3plGjRpx22mkly1999dUsXry4ZDhos2bNGDduXKJ+XSDE6SXNbD2wi3POmW8LrXfOlTtWsqrTSyZ1+r8lS/yIhzff9Hf9jhjhb/wSEQDmzp1L27Ztww4ja5Sxv9NyesnlQCfgE+AUYEGyNpSU6f+2bNlW6dM5eOABuOkmVfoUkYwSZhL4I/CYmeUA+cBVIcZSOZ9+6sf8z5mjSp8iktFCSwLOuanAEWFtv0pWrPBj/keP9gf9d95RoTcRyWiqHRQPVfoUkRoqK8pGVMvXX/uunxkzoGtXX+mzTZuwoxIRSQi1BMqydq0v83zssdsqfU6YoAQgIjWKkkAs5+DFF/3BftQoX/Zh3jy4+GLQxBgiGal27dp06NCBgw8+mPbt2/Pwww9TXFxc7jqLFi3ilVdeSVGE4VF3ULTvv/dn/1OnqtKnSEiSUfW3QYMGJXf5rly5kksuuYT169czaNCgMteJJIFLLrmkWttOd2oJAGzYAP37+wP+3Lnw3HOq9CkSglRU/W3WrBmjRo3iiSeewDnHokWL6NixI4cffjiHH354SVmH2267jSlTptChQwceeeSRMpfLdNndEnAO/vEPf8fvsmXwxz/6Sp9Nm4YdmUhWSlXV33322Yfi4mJWrlxJs2bNmDhxIvXr12fBggV0796d6dOnM2TIEIYNG8Z7770HQF5eXqnLZbrsTQK5udC3L3z4oT/jf+MNfxFYREKTyqq/kZI5hYWF9O3bl2+++YbatWvzn//8p9Tl410u02RfEti8GYYOhSFDfImHxx7zQ0Bzsm9XiKSbSNXf6ESQjKq/P/zwA7Vr16ZZs2YMGjSI5s2b8+2331JcXEz9+qUnm0ceeSSu5TJNdl0T+Oc/fZnnQYPg//7Pj/rp108JQCRNRKr+1svxh6ZkVP1dtWoV11xzDX379sXMWL9+PXvuuSe1atVi9OjRJTN7NW7cmA0bNpSsV9ZymS47jn7OwUUX+cldDjwQ/vUvOOWUsKMSkRgN6+XQtW1zJg/onNCqv5s3b6ZDhw4UFhaSk5NDjx49uOmmmwC49tprOf/883n99dfp3LlzyUQz7dq1Iycnh/bt29OzZ88yl8t0oZWSroqqlpIG4K67oFEjX+kzjtmBRCQxVEo6tTKplHRqBRM6i4jINtl1TUBERLajJCAiSZdJ3c6ZrCr7WUlARJKqfv36rFmzRokgyZxzrFmzptJDV7PnmoCIhKJly5YsXbqUVatWhR1KjVe/fn1atmxZqXWUBEQkqerUqcPee+8ddhhSBnUHiYhkMSUBEZEspiQgIpLFMuqOYTNbBfwUdhyVsBuwOuwgQqZ94Gk/aB9EhLEfVjvnupX2QkYlgUxjZtOdc0eGHUeYtA887Qftg4h02w/qDhIRyWJKAiIiWUxJILlGhR1AGtA+8LQftA8i0mo/6JqAiEgWU0tARCSLKQmIiGQxJYEkM7OHzGyemc0ys7fMbJewY0o1M7vQzGabWbGZpc3QuFQws25mNt/Mcs3strDjCYOZPW9mK83s+7BjCYuZ7WVmH5vZ3OB/4YawY4pQEki+icAhzrl2wH+A20OOJwzfA/8HTA47kFQys9rACOAM4CCgu5kdFG5UoXgRKPVGpSxSBNzsnGsLHAtcly6fBSWBJHPOTXDOFQUPvwQqV+e1BnDOzXXOzQ87jhAcDeQ6535wzm0BxgLnhBxTyjnnJgP/CzuOMDnn/uucmxn8vAGYC7QINypPSSC1rgQ+CDsISZkWwJKox0tJk398CY+ZtQYOA74KNxJP8wkkgJlNAvYo5aU7nHPjgmXuwDcJx6QytlSJZx9kISvlOY3JzmJm9ivgH8CNzrlfwo4HlAQSwjnXtbzXzexy4HdAF1dDb8yoaB9kqaXAXlGPWwLLQ4pFQmZmdfAJYIxz7s2w44lQd1CSmVk34FbgbOdcXtjxSEpNA/Y3s73NrC5wMfBOyDFJCMzMgOeAuc65h8OOJ5qSQPI9ATQGJprZN2b2VNgBpZqZnWdmS4HjgPfN7MOwY0qFYEBAX+BD/IXA15xzs8ONKvXM7FXgC6CNmS01s15hxxSCE4AewCnBceAbMzsz7KBAZSNERLKaWgIiIllMSUBEJIspCYiIZDElARGRLKYkICKSxZQERESymJKAiEgWUxIQEcliSgIiVWRmu5vZf83s7qjn2plZvpldEGZsIvHSHcMi1WBmpwPvAp2Ab4DpwNfOuStCDUwkTkoCItVkZo8CZwOfAh2BDs65jeFGJRIfJQGRajKzesC3wP7A8c65tJgsRCQeuiYgUn2t8fMGOGCfcEMRqRy1BESqIZgo5AtgAX66wHuBds65xWHGJRIvJQGRajCzIcAlQDtgPX4O6QZAZ+dccZixicRD3UEiVWRmnYCbgcucc+uCqUN7Am3xs8mJpD21BEREsphaAiIiWUxJQEQkiykJiIhkMSUBEZEspiQgIpLFlARERLKYkoCISBZTEhARyWL/D+PoBOITC9ARAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "# generate data\n",
    "np.random.seed(123)\n",
    "N = 20\n",
    "beta0 = -4\n",
    "beta1 = 2\n",
    "x = np.random.randn(N)\n",
    "e = np.random.randn(N)\n",
    "y = beta0 + beta1*x + e\n",
    "true_x = np.linspace(min(x), max(x), 100)\n",
    "true_y = beta0 + beta1*true_x\n",
    "\n",
    "# plot\n",
    "fig, ax = plt.subplots()\n",
    "sns.scatterplot(x, y, s = 40, label = 'Data')\n",
    "sns.lineplot(true_x, true_y, color = 'red', label = 'True Model')\n",
    "ax.set_xlabel('x', fontsize = 14)\n",
    "ax.set_title(f\"Linear Regression for y = {beta0} + {beta1}x\", fontsize = 16)\n",
    "ax.set_ylabel('y', fontsize=14, rotation=0, labelpad=10)\n",
    "ax.legend(loc = 4)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<u>Parameter Estimation</u>**\n",
    "\n",
    "The above covers the entire structure we assume our data follows in linear regression. The machine learning task is then to estimate the parameters, $\\bbeta$. These estimates are represented by $\\hat{\\beta}_0, \\dots, \\hat{\\beta}_p$ or $\\bbetahat$. The estimates give us *fitted values* for our outcome variable, represented by $\\hat{y}_n$. \n",
    "\n",
    "\n",
    "This task can be accomplished in two ways which, though slightly different conceptually, are identical mathematically. The first approach is through the lens of *minimizing loss*. A common practice in machine learning is to choose a loss function that defines how well a model with a given set of parameter estimates approximates the observed data. The most common loss function for linear regression is squared error loss. This says the *loss* of our model is proportional to the sum of squared differences between the true $y_n$ values and the fitted values $\\hat{y}_n$. We then *fit* the model by finding the estimates $\\bbetahat$ that minimize this loss function. This approach is covered in the first sub-section.\n",
    "\n",
    "\n",
    "The second approach is through the lens of *maximizing likelihood*. Another common practice in machine learning is to model the outcome $y_n$ as a random variable indexed by one or more parameters and then find the parameters that maximize its likelihood. The most common model for $y_n$ in linear regression is a Normal random variable with mean $\\bbeta^\\top \\bx_n$. That is, we assume \n",
    "\n",
    "$$\n",
    "y_n|\\bx_n \\sim \\mathcal{N}(\\bbeta^\\top \\bx_n, \\sigma^2),\n",
    "$$\n",
    "\n",
    "and we find the values of $\\bbetahat$ to maximize the likelihood. This approach is covered in the second sub-section. \n",
    "\n",
    "\n",
    "Once we've estimated $\\bbeta$, our model is *fit* and we can make predictions. The below graph is the same as the one above but includes our estimated line-of-best-fit, obtained by calculating $\\hat{\\beta}_0$ and $\\hat{\\beta}_1$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEcCAYAAAAhoQi5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3gU5fbA8e9JQhICoRN6EwFBqSJNWigBQZqKCIpYsIBcAQWkeBVEFBQEvYCKXjvKVa6AP70mIbQAUkITQUBAOgoECSSE9Pf3x2xiCOnJZpLs+TwPT9jdKWdnd+fMvPPOecUYg1JKKdfkZncASiml7KNJQCmlXJgmAaWUcmGaBJRSyoVpElBKKRemSUAppVyYJoFcEJFHRMSIyM2ZTPOJiBwvwLDyhYjUdby35H+JIvKniCwVkVp2x1cQROS4iHxiw3rbi8g2Ebnq2PYtCjqGwkZEOohIkmN7eBTgehuKyNsisldEokTkDxH5TkSaF1QMBaXANqoLmgm8bXcQefA68B3gCbQDXgYai0hbY0y8rZE53yDgig3r/TdwDegHRAO/2RBDoSEiJYD3gXNA1QJefQDgD3wK7ALKAZOAbSJypzFmZwHH4zSaBJzEGHPU7hgy4vhxJZjM7xT83Riz1fH/UMc8rwK3A1szni1/iYg7IMaYhIJapzFmd0GtK5mIuAGNgFnGmLX5sDwBShhj4vIcnH0mAgJ8BEzNzQJEpC5wDPA3xqzPwazLgEWpfyMishY4DowFHs5NPIWRNgc5SdrmoFTNLE+JyCuO08sIEfk/EamZzvxPiMjPIhIjIuEi8m8RqZBmmjEiskVE/nIsa6uI9E0zTfJ6R4vIGyJyFojFOrLJiV2Ov7XTLN9HROaIyDERiXP8nebYqaWerpWIbHS8n1MiMlVEZoiISTOdEZFZIjJZRI4BcUBTx2uVRORdETkjIrEiclBEnkwzf1UR+VREzjqm+UNEvhcRP8frHiIyU0SOptq2m0SkY6pl3NAcJCJtRCTE0TRwVUTWiEibNNN8IiKnRaSl471Gi8hhEXk6sw0rIo8AiVi/x386tsHxVK8/lOa78LmIVEuzjOMi8oWIPCYiBx3b7brvQqppfxGRFek839Wx7l6ZxVsQRKQ+MA0YDRT4macxJjztQZIx5jLW2VmN5OdEZKRjmw1M9Zy7iIQ6vmO+BRd1Lhlj9F8O/wGPAAa4OZNpPgGOp3pc1zHPceBL4C5gBBAObEgz72ysL/48rNPSR4EzwDbAPdV0c4HHge5AL2ChYx13pbPeM8BK4G5gAFAyg7iTpx+Z5vlRjudvT/WcB7ARuAiMc8QxDYgB5qWarhJwCdgP3A8MBDYAJ6yv4HXrSY51I3Av0BuoApQBDgEngSeAHsCbWDvPf6SafzXWD/VBoDMwGHgPqOt4fRoQhXU01wWr6WUG0D/VMo4Dn6R63AyrmWYncJ8jrjDHc83TfOZXgAPAU0BPx2dtsI5EM/quVAbudEz3IVbzW0vHa086nl8G9AFGAucd77F0mpjPAPuAoY7Pon4G60vesVZP8/xXwO9YZ16Zff89svHPPbNlZOM3Fgx86vj/dMc28MjFcuo65u2aD7/7CsBVYGGa57/G+h3XSBVvPNA2r+ssiH+2B1AU/5G3JJB2hz/B8Xz1VNMlAi+lmS55JzEwg/W5OX58wcCqdNa7K6sfd5rpn3QszwfoBpwGlqeZdrhj2s5pnp+GdSTq53j8muNxzVTTlMRq6zVp5jXAWdIkKeCfWMmlQZrnP3D8AD0cj6OAZzN5f98D32axDY5zfRJYDkQA5VI9Vwb4K/WyHJ/5dTt8wMsR35Is1unhmHd6qufcHdtoXZppOzqmfTZNzNFA1Wx8xr5YyeqfqZ6rhHWGODmb34+s/h3PKo5M1vGQY9smf3+mk80kkOp3kPyvvmPe7mmed8tFXEsd2/jmNM+XwzqgWYd1YJEATMnt+y/of9ocVPB+SPP4F8ff5GaWnlhf5KWOpgsPsXpFbMP64XZOnlFEbnc0dZzD+uLFO+ZvlM56VxrHNzab3ncs7yqwBmtn9FCaaXpjffl/ShNrMFAC64gWx98txpjTyTMaY65x47ZIFuh4Pe26tgHH0qwrCKgINHFMFwZMFJGxItLU0TaeWhjQx9Hk1FFEPLOxLToD3xtjIlLFfwXrwnmXNNNGG2PWpZouFjhMmma0bGoE+GHtfFIYYzZhbfe0695qjPkzq4UaYyKBL4CRqZrtHsVqf/84i9nPAndk41+/zBYiIm6pP0exrv3gaPKcB0w1xpzP6r2k4yWs723yvyOO50PSPP9SThYqIlOAYcAYY8yR1K85vhfDgE5Y38eNwJxcxG4LvTBc8P5K8zjW8dfb8dfP8fcI6asIIFZ3zTXAr8A/sJpJErB6JTVOZ74/chjnq8AqrDOBwcAYYDHwWKpp/IA6ZNxmW9HxtxpWM0Va5zKYL71Y/YCbs7GuIVg9mSYBC4A/ROQ94FVjTBLWWUkMVkKbCkSJyHJgojEmPINlV8ggpj+B8mmeu5TOdLH8/fnmRPI1oIzWXSHNczn5jBdjNfH1EZEfsM78VhhjMvpMADDGxInInmwsP6sDjpewPqdkG4CuWN+7c8DXIpJ83Sp525UVkRhjzNVMlrsE62wvWTWsZP00VnNesrNZxJfCcU3nNeBFY8xHGUy2Fau5sgnwtuO7ViRoEih8Ljr+BpD+DiX59d5AWeD+1EfYIuKTwXJzWjP8hDFmh+P/oY4LXI+KyHvGmO2pYjmG1c6fnuOOv3/wd3JLrUoOYr2I1RY+NoN5DgE4jh6fAZ4RkUZY111mABeAd43VvXUOMEdEqmJdI3kLK9kNyWDZf5F+F8Wq3JjU81PysjNa9440z2X7MzbG7BORjVjXLmKwEuxTWc0nf/e2ycoJrKajjKTdWUc6/jbB6ghw8YY5rGa1VVjXlNJljDlLqh28I16AQ6m+z9kmIsOxEuY8Y8ysTCZ9GWgA7AXmi8g6Y11ILvQ0CRQ+q4EkoLYxZnUm0yXv7FOOjEWkIda1g9PpzpE3k7F29i/zd6+TQKyLpFHGmIOZzLsVmCAiNZMTloiUJIPeKxkIxHHGk91mAmPMIWCq40jutnRe/xP4UET6pPd6KhuAviLi62hKwZEU+wHrc/AecuoQ1lHxA1j3EOBYdwesM7B5eVz+YqxmofLAbyZ7XVOTm4OyEpvZi2l31qmM48aea49gJfMeZHz2mO9EZBBW89iHxpgJmUzXCeuscgrwH+Bn4F2sJqJCT5NA3vQWkbRtsJez2HlnyhhzVETmAAsdR7IbsI7UamG193/oaHMOwWr++UxE5mGd9s7AahbK92s9xpg/RWQR1s78dmPdLLMUqy15jSOGn7FuLqsP9Me6iB2NdaQ9CggSkRlYO4jnHH+ze/Q6H+tIfaOIzMfaQZYCbgE6GWMGiEhZrO2yFDiIlSAHYO3kggFEZJUjzl1YZ1otsc6q3s9k3TOxzhjWOD4bA7yAlYhfyWb8OWaMSRSRl4D3ReQLrB12DWAW1nWGrNrvs/JfrCazO4HnsxlTHDeegeQbY8wNTU0i0tXx3w2mgO4XEZHOWL2l9gKfiEi7VC/HGse9JCJSHuv7tg6Ya4wxYnVb/lpEgowxnxZEvHmhSSBv/pXOc/vJ/KgyS8aYqSJyAEezBtZO5xTWNYDDjmn2i8iDWDuh74CjWEfrvbHaVp1hNlaTwUvAAGNMvKNP+WSsNuV6WBeSj2Jd9I1zxBouIt2Bd4DPsE7138PqkZKtm26MMZcdR8AvYe2Aa2D12DmEtTMDK1nuwupCWgfrjOoQ8KAxZpVjmlCsaxzPYO3ETwJvYO1YM1r3XseOaBbWHaSCdXbTxRjzc3bizy1jzBIRica6cWoVVu+n/wGTjDFReVx2vCMpjsB6X+pv3bB6drUENqd5LXVT1xKsnm4PJ3e8MMZ8IyL/xjqQ25z2QnJhIznrMKJU/nD0BtkFhBtjutsdjyty9K46Amw0xgy3Ox5lDz0TUAVCRGZi7XBOYPXkGYl1E1YfO+NyRSJSButsdRhWM2Nery2oIkyTgCooBqspp7rj/3uxrhn8aGtUrqkVVhv2eWBseu3wynVoc5BSSrkwvWNYKaVcWJFqDurdu7cJDAy0OwyllCpq0pZPSVGkzgTCwzO6q18ppVRuFKkkoJRSKn9pElBKKRemSUAppVyYJgGllHJhmgSUUsqFaRJQSikXVqTuE1BKKVcRHZtAZGwCl6LjKO/jia+XBz5e+b/L1iSglFKFTHRsAiEHzjHhm73EJSbh6e7G3MHN6NG4Sr4nAm0OUkqpQiYyNiElAQDEJSYxcfleImPzf0wdTQJKKVXIXIqOS0kAyWITkoiIjs9gjtzTJKCUUoVMeR9PPN2v3z17ebhRzqdEvq9Lk4BSShUyvl4ezB3cDC8Paxft5eHGm/c1x9cJF4aL1HgCrVu3Njt2OG2Ma6WUKjSSewdFRMdTzqdEXnsHZVhFVHsHKaVUIeTj2OlXKePt1PVoc5BSSrkwTQJKKeXCNAkopZQL0ySglFIuTC8MK6VUIVUQ9YM0CSilVCFUUPWDtDlIKaUKoeT6QX4Xz3LruaNOqx+kSUAppQqhiEuRPLnxS0L+PZpZQYvAGKfUD9IkoJRShU1ICPV73MmEjV8QUr8NTw+cCiJOqR+k1wSUUqqwOHsWnnsO/vMf3OrfzJZFX/D82QrEJiQ5rX6QJgGllLJbQgL861/w8ssQFwczZuA2aRLNxYPQ/KsflC5NAkopZafNm2H0aNi7F+66y0oG9esD4ANOrx+k1wSUUsoOFy7A449Dx45w6RJ8+y388ENKAigotiYBEektIodE5IiITLYzFqWUKhBJSbBkCTRqBJ99Bi+8AAcOwKBBIBlWfHYa25KAiLgDi4C7gCbAUBFpYlc8SinldLt2QYcO8NRT0KwZ/PwzzJ4NpUrZFpKdZwJtgCPGmN+NMXHAMmCAjfEopZRzXL4Mzz4Ld9wBx47B55/DunXQxP7jXjuTQA3gVKrHpx3PXUdEnhSRHSKy48KFCwUWnFJK5ZkxsHSp1fSzaBGMGgWHDsFDD2W76ccYQ1xUnNNCtDMJpLcFbhjr0hizxBjT2hjTunLlygUQllJK5YMDB6B7d2uHX7s2bN8OCxdCuXLZXsSZsDN8dOdHfPf4d04L084kcBqolepxTeCsTbEopVT+uHoVJk+22vz37IH33oMtW+D227O9iMizkawcsZIP23zIpd8vUb+383oM2XmfQBjQQETqAWeAB4BhNsajlFK5ZwysWgVjx8LJk/DIIzBnDvj5ZXsR8dfi2fLWFja9vomk+CTunHwnnaZ2wsvXy2lh25YEjDEJIjIGCALcgY+MMfvtikcppXLt99+tC78//AC33QahodCpU7ZnN8bw6/JfWT1xNZdPXKbxPY3p8UYPKtSv4MSgLbbeMWyM+R/wPztjUEqpXIuNhTfegNdeAw8PmDvXSgYlsl/k7Y/dfxA4NpCTG09SpVkVBqwdQD3/ek4M+npaNkIppXJj9Wp45hk4fBjuuw/mz4eaNbM9e9SfUax9cS27P9qNTyUf7n7/blo+3hI394K9VKtJQCnlVAUxRGKBOnPGqvT59ddw880QGAi9emV79oTYBLa9vY3QV0NJiEmg/XPt6fzPzniXdV59oMwU4U9CKVXYFdQQiQUiPv7vSp/x8TBjBkyaBN7Z23kbYzi48iCrJ6zm0u+XaNivIQHzAqjYoKKTA89cEfsUlFJFSfIQiXGJSQApQySGTvIvWklg82brRq9ffoE+feCdd3JU6O3c3nMEjQ/i2NpjVG5SmYeCHqJ+QMEWistIEfoUlFJFzaXouJQEkCx5iERnlkfONxcuWAXePv4YatWyKn0OHJjtu32vXrjKun+uY9cHu/Au581dC++i9VOtcfMoPAWcNQkopZymvI8nnu5u1yUCZwyRmO+SkuCDD2DKFIiMtJp9Xnop24XeEuMS2b5oOxtmbCAuKo47xtxB15e7UrJCSScHnnOaBJRSTuPr5cHcwc2YuHyvU4dIzFe7dllNP9u3Q9euVs2fbBZ6M8Zw+IfDBD8fzMXfLnJz75sJeCuAyo0Lb8mbQvxJKKWKOh8vD3o0rkLoJH+nDpGYLyIi4J//hMWLoXJlq9Lngw9mu+nnwq8XCHouiKNBR6nYqCLDfhhGgz4NnBx03hXCT0IpVZz4OHb6hfYagDHw5Zfw/PPWNYDRo2HmzGwXeou+GM366evZ8e4OPEt70mt+L+545g7cS7g7OfD8oUlAKeW6Dhywdvrr10ObNlbZh2wWekuMT2THeztY//J6Yi/HcvtTt+P/ij8+lXycG3M+0ySglHI9V6/Cq69aZR58fa1Kn088AW7Z67VzJOgIQeODCD8QTr3u9eg1vxdVmlZxctDOoUlAKeU60lb6fPRRq9JnNscqCT8UTvDzwRz+4TDl65dnyMohNOrfCLFhbOD8oklAKeUa0lb63LgROnbM1qwxETFseGUD2/+1HY+SHvR4owdtn22LR2G8wJ1DRf8dKKVUZtJW+pw3D/7xj2xV+kxKTGLXB7tY9891RF+MpuXjLen2ajdKVyldAIEXDE0CSqniK3Wlz8GD4a23sl3p89jaYwSOC+T8L+ep07kOvRb0olrLak4OuOBpElBKFT95qPT519G/WD1hNQdXHqRc3XIM/mYwje9tXKTb/TOjSUApVXykrvSZkJCjSp+xV2IJnRXKtgXbcCvhRrdZ3Wj/XHs8vIv3brJ4vzullOvYtMnq859c6fNf/4KbbspytqTEJPZ8soe109Zy9dxVGg67jVsmtqf6TRWIk+K/kyzu708pVdylrfS5YgUMGJCtcg8nQk8QOC6QP3f/SfV2Nakxrwcv//oHcct2F+2xD3Kg8NQzVUqpnEhKgvffh0aNrDo/kyZZdwBno9RzxPEIvhn8DZ90+YTo8Gju/epe+gUOsxJAmrEPImMTCuLd2Kb4pjelVPG1c6fV9JPDSp9xUXFsfH0jW+Ztwc3dja4zutJhQgdK+JTg4J9XivbYB7mkSUApVXSkrfT5xRcwbFiWR/4mybD3i72ETA4h6o8omj7YlB6ze1CmZpmUaYrs2Ad5pElAKVX4GQNLl8KECTmu9Hnqp1MEjgvkbNhZarSpwZBvh1Cz3Y33ChTJsQ/yQfF+d0qpou/XX0l8ehTuG0O51vJ2or9eQcm2d2R5sfbyqcusmbyGX778Bd/qvgz8bCDNHmyGuKV/1lCkxj7IR8X73Smliq6rV2HmTMy8eST6lGL6Xf/gi6Y9KREUwdwy5zLstRMfHc/mNzezec5mMNDpxU50fKEjnqU9s1xloR/7wAk0CSilCpfkSp/PPgunThHz0MN0rXgX57x9gb977YRO8r8uCRhj2LdsHyEvhHDl1BVuvf9WerzRg3J1sjc4jKvSJKCUKjxSV/ps2hS++ooT9ZtybsHG6yZL22vnTNgZgsYFceqnU1RrVY17lt5DnU517HgHRY4mAaWU/TKp9Fn+SkyGvXYiz0ayZuoafv70Z0pVKUX/f/en+YjmuLnrLVDZpUlAKWWv4GAYM8aq9Hn//Valzxo1Ul5Or9fOnH63sn/BNra+sZmk+CTufOFOOk3thFcZLxvfSNGkSUApZY/UlT4bNICgIAgIuGGy1L12Ll2N46/g39l2/7dcOXmZWwbeQs+5PalQv4INb6B4sCUJiMibQD8gDjgKPGqMibAjFqVUAUtb6XPmTJg4EbwyPor38fLg8q8X2DI2kJMbT1KlWRUGfjyAet3qFWDgxZNdZwKrgSnGmAQRmQNMAV6wKRalVEFJXemzb194550sK31GnYti7bS17P5oNz4Vfej7Xl9ajWyl7f75xJYkYIwJTvVwK3CfHXEopQrIhQtWgbdPPoHatWHlSujfP9NyDwmxCWx7exuhr4aScC2BduPb0eWfXfAu5zp9+AtCYbgm8BjwH7uDUEo5QWIifPghTJkCkZEweTK8+CKUKpXhLMYYDq06RPCEYC4dvUTDfg0JmBtAxYYVCzBw1+G0JCAiIUDVdF6aZoxZ5ZhmGpAALM1kOU8CTwLUrl3bCZEqpZwiF5U+z/1yjqBxQRxbe4zKTSrzUNBD1A+oXzDxuigxxtizYpERwNNAd2NMdHbmad26tdmxY4dzA1NK5U1EhHW0/+67VqXPt96CoUMzbfq5euEq615ax64lu/Au503XGV1p/XRr3Dy03T+fZLjx7eod1BvrQnCX7CYApVQhl1zp8/nnITw8W5U+E+MS2b5oOxtmbCAuKo47nrmDrtO7UrJCyQIM3LXZdU1gIeAFrBbr6GCrMeZpm2JRSuXVr7/CM8/A+vXQpg38+CO0apXh5MYYDv/vMMHPB3Px0EXq96pPr7d6UblJ5YKLWQH29Q662Y71KqXymaPSJ/Pmga+vNdzjyJHglnEzzoVfLxD0XBBHg45SsWFFhn4/lAZ9GiDZGBNY5b/C0DtIKVXUpKn0yWOPwezZ1jWADFz76xrrp68nbHEYnqU9CXgrgDbPtMHd070AA1dpaRJQSuVMOpU+ufPODCdPjE9k5/s7Wf/yemIiYmj1ZCv8X/GnVOWMu4mqgqNJQCmVPWkrfb71llXp0yPj3ciRoCMEjQ8i/EA49brVo9f8XlRpVqUAg1ZZ0SSglMpa6kqfQ4ZY1wBSVfpM6+JvFwl+Ppjfvv+N8vXLM2TFEBoNaKTt/oWQJgGlVMbOnIHx4+GbbzKt9JksJiKGDa9sYPu/tuNR0oMec3rQdmxbPIr5OL1FmX4ySqkbpa30+corVu2fDCp9JiUmsevDXax7cR3RF6Np+XhLur3ajdJVShdw4CqnNAkopa6XutJnnz5WMsik0uextccIHBfI+V/OU6dzHXot6EW1ltUKMGCVF5oElFKWHFb6/OvoX6yeuJqDKw5Srm45Bn8zmMb3Ni507f7RsQlExiZwKTqO8j6e+Hp5XDdAvavTLaGUq0tMhA8+gKlTs1XpMzYylo2zNrJ1/lbcSrjRbVY32j/XHg/vwrc7iY5NIOTAOSZ8s5e4xCQ83d2YO7gZPRpX0UTgoFtBKVe2cyeMGgVhYVlW+kxKTGLPJ3tYO20tV89dpfmI5nR/rTu+1X0zXYWdR+KRsQkpCQAgLjGJicv3EjrJX5OAg24FpYq4XO1kkyt9Ll4Mfn7w+efw4IMZNv2c2HiCwLGB/Ln7T2q2r8nQ74ZSo03GXURTx2bnkfil6LiUBJAsNiGJiOh4qpTRwWlAk4BSRVqOd7JpK32OGWP1/Mmg0mfE8QhWT1rNr9/8SpmaZbjny3u47YHbst3ub/eReHkfTzzd3a5LBF4ebpTzKeH0dRcVWqxbqSIso51sZGzCjRPv3w/+/jB8ONStazUBvfNOugkgLiqOtS+uZeEtC/nt+9/oMr0LYw6NoenQpjm68JvZkXhB8PXyYO7gZng5xiXw8nDjzfua46tNQSl0SyhVhGWruSMqyqr0+dZbVqXPJUvg8cfTrfRpkgx7v9jLmilriDwbSdNhTek+uztla5XNVXx2H4n7eHnQo3EVQif5ExEdTzmfEto7KA3dEkoVYZnuZI2xunmOHZutSp+ntpwiaFwQZ7afofod1Rm8fDC12tfKU3zJR+ITl+8lNiHJliNxH8dOX68BpM+24SVzQ4eXVOp6ydcE0u5ke3pGUnLCeGtwl6ZNraEeM6j0efnUZdZMXsMvX/5C6Wql6TG7B80eaoa45U9//+QL13okbqsMP0xNAkoVcdftZN0SKb9wAZ5vzIYSJayLvhlU+oyPjmfzm5vZPGczGGg/oT0dX+iIZ2lPG96FcrLCNcawUir/pDR3bA21hng8cgTuv9+6BpBOpU9jDPuW7SPkhRCunLpCk8FN6PlGT8rVzXgsYFV8aRJQqqg7fRqee+7vSp/BwdCzZ7qTnt1xlsCxgZz66RRVW1blni/uoU7nOgUcsCpMNAkoVVSlrfQ5cyZMnJhupc/Is5GsmbqGnz/9mVJ+pej3YT9aPNICN3ftJe7qNAkoVRSlrvTZt6+VDOrVu2GyhJgEtry1hY2vbSQpPokOkzrQeVpnvMqkXxJauR5NAkoVJdms9GmM4cB/D7B64moijkdwy8Bb6Dm3JxXqV7AnblVoZZkERORhYD5Q3RgTm+r5pYCvMaa/E+NTSoFV6fPDD2HKFOvmr0wqff6x+w+CxgVxIvQEVZpV4eE1D1Ov241nCUpB9s4EvgHeBgYAXwOISFlgEDDUeaEppQCr0ufo0bB9u1Xpc/FiaNz4hsmizkWx9sW17P73bnwq+tD3vb60GtlK2/1VprJMAsaYa46j/sdwJAFgGHAF+MGJsSnl2pIrfb77rnWX79KlMHToDU0/CbEJbHtnG6EzQ0m4lkC78e3o8s8ueJfTO2RV1rJ7TeADYJeI1DTGnMZKCJ8aY9KpUqWUyhNj4IsvYMIEq9Ln6NFWz580hd6MMRxadYjgCcFcOnqJhnc3JGBeABUbVrQpcFUUZSsJGGN+FpFdwCMishJoDTzk1MiUckX791s7/dBQaNvWKvvQqtUNk5375RxB44M4tuYYlZtU5qGgh6gfUN+GgFVRl5PeQR8Ak4BKwGZjzCHnhKSUC0pd6bNMmQwrfUaHR7PupXXsfH8nXmW9uOtfd9H66da4eWi7v8qdnCSBr4C3gFHA084JRykXYwysWGFV+jx92trxz54NlSpdN1lifCJhi8LYMGMDsZGx3PHMHXSd3pWSFUraFLgqLrKdBIwxkSLyNTCYvy8QK6Vy6+hRePZZ+N//oFkzWLYs3Uqfh/93mKDngrh46CL1A+rTa34vKjdJvxy0UjmV05vFqgHLjDFXnRGMUi4hJgbeeANee82q9Dl/vjXMY5pKnxcOXCD4uWCOBB6hYsOKDP1+KA36NMjRyF5KZSVbSUBEKgA9gACgeX4GICITgDeBysaY8PxctlKp5WpA9vwWFGTt8I8cgSFDYN68Gyp9XvvrGutnrCdsURiepT0JmBdAmzFtcPd0L9hYlUvI7i9gF1ABmGqM2ZdfKxeRWkBP4GR+LVOp9OR4QPb8dvo0jB8Py5dnWOkzKSGJHe/tYP3L64mJiKHVE63wn+lPqco33hWsVH7JbgNfuCcAACAASURBVBfRuk5a/3ysHkernLR8pYCMB2QPneTv3CQQH28N5j59ulXp89VXrf7/aSp9Hg0+StD4IC78eoF63erRa34vqjSr4ry4lHKwrYCciPQHzjjuQchsuieBJwFq165dQNGp4iZbA7Lnt02bYNQo2Lcvw0qfF3+7SPDzwfz2/W+Uv6k8Q1YModGARtrurwqMU5OAiIQAVdN5aRowFesaQ6aMMUuAJWANL5mvASqXkemA7Pnt/Hl44YVMK33GRMSwYeYGtv9rOx7eHvSY04O2Y9vioWPvqgLm1G+cMaZHes+LSFOgHpB8FlATqyxFG2PMn86MSbkmXy8P5g5udsOA7L75udNNTIQPPrAqfV69mm6lz6TEJHZ9uIt1L64j+mI0LR9rSbdZ3ShdpXT+xaFUDhSKgeZF5DjQOqveQTrQvMqL6wZk9ymRv72Ddu60mn7CwsDfHxYtuqHS57F1xwgaF8S5veeo3ak2vd/uTbWW1fJn/UplTgeaVyplQPb8vAaQXOlz8WLw87MKvw0bdl3Tz6XfLxE8IZiDKw5Stk5Z7vv6Pprc10Tb/VWhUCiSgBN7HynlHGkrfY4ZY9X+KVs2ZZLYyFg2ztrI1vlbcSvhhv+r/rR/rj0lSub8OkShuMdBFUv6LVIqp7Ko9GmSDHs+2cOaqWu4eu4qzR9uTvfXu+Nb3TdXq7P9HgdVrGnpQaWyKyrK6vXTooXV7XPJEvjpp+sSwMlNJ/ngjg/47vHvKH9TeUZuG8nATwfmOgFAxvc4RMbqcB4q7/QwQqmsZKPSZ8SJCEJeCGH/f/ZTpmYZ7ll6D7cNvS1f2v1tucdBuQxNAkpl5sgRq9Lnjz9alT7/8x/o0CHl5bioODbN2cSWuVtAoMvLXegwsQOepTzzLYQCvcdBuRxNAkqlJyYG5syB119Pt9KnSTLsXbqXNZPXEHk2kqbDmtJ9dnfK1iqbxYJzrkDucVAuq1DcJ5Bdep+AKhCBgdYO/+jRdCt9nt56msCxgZzZfobqd1Sn99u9qdW+llNDcuo9DsoV6H0CSmUpdaXPhg1h9Wro8fdN71dOXyHkhRB++fIXSlcrzcBPB9LsoWaIm/P7+zvlHgel0CSglFXp8+23rUqfiYlWf/+JE1MqfcZHx/PT3J/YNHsTJsnQaVonOk7uiGfp/Gv3V8oumgSUa9u40Sr3sH8/3H23VfbZUenTGMP+/+xn9aTVXDl1hSaDm9DzjZ6Uq1vO5qCVyj+aBJRrOn8eJk2CTz+FOnVg1Sqr0qfD2R1nCRwXyKnNp6jasir3fHEPdTrXsTFgpZxDk4ByLWkrfU6ZAtOmpVT6jPwjkrVT17Lnkz2U8itFvw/60eLRFri5632VqnjSJKBcRyaVPhNiEtgyfwubXttEYlwiHSZ1oPO0zniV8cpioUoVbZoEVPEXEQHTpmHefZekyn6cW/gB7g89iK93CUoaw4H/HmD1xNVEHI/gloG30PPNnlS4uYLdUStVIDQJqOIrVaVPEx7OsSGPcF/1Pvx1qiSeb6znlabVubpoF6c3ncSvqR8Pr3mYet3qZb1cpYoRTQKqeEpT6fOv5avoHXiJuMQkvK/G0yr0DAdf207JiiXp+25fWo1shZuHtvsr16PfelW8REVZvX7SVPq80KAJCbEJ3LbtT+5d8gsN9l1kf+sq9Nn8KK2fbq0JQLksPRNQxYMx8O23MG7cDZU+jTFErDvBoI/2U+ZSLCfrlyWsWy1i/Xzwq5b7Es9KFQeaBFTRd/SoVesnMBCaN7+u0ue5X84RND6IY2uO4VevHD/2qsPxOmW0CJtSDvoLUEVX6kqfnp6wYAE88wx4eBAdHs26l9ax8/2deJX14q5/3UXjR5rzSJLRImxKpaK/AFU0pa70+cADVqXP6tVJjE8k7O2tbJi+gdjIWFqPbk3X6V3xqegDgC9oETalUtEkoPKsQAdBP3XKqvT53//eUOnz8P8OE/RcEBcPXaR+QH0C3grA71Y/58ShVDGhSUDlSYENgp620uesWfD88+DlRfjBcILGB3Ek8AgVGlRg6P8NpUHfBvkytKNSxZ0mAZUnGQ2CHjrJP/+SQAaVPq/9dY31k34kbFEYnqU9CZgXQJsxbXD3dM+f9SrlAjQJqDxx6iDoaSp9xi7/loied3HxSgx/zN3Mjtc3ExsRQ6snWuE/059SlUvlbX1KuSBNAipPnDIIemKidZPX1KkplT6jJ7xAyMko5j61klarT1D+Ygyl21RnyKI+1G1dI+tlKqXSpbdJqjxJHgTdy3HHbZ773+/cCe3bWyUfWrWCvXvhtdc4dfwq/3f/crovO4R7omHNoPos6VGDkg0r5uO7Ucr16JmAyhMfLw96NK5C6CT/vPW/v3QJXnwR3n0XqlSBpUth6FBiLsey4fkgtv9rO34CYV1r8OvtVUjycINEkz/NTkq5ME0CKs/yNAh6qkqfhIfDP/4Br7xCUmlfdi3ZyboX1xF9MZpbHmrGrPJuXCn591c2z81OSilNAsVRgfbbz4s0lT4JDISWLTm+/jiB477i3M/nqN2pNr0X9KbsrZUpeeAcE5fvJTYhScs+KJVP9BdUzBRYv/28iIqCV16B+fOhTBlruMfHHuPS8csE3/MfDq44SNk6Zbnv6/tocl+TlP7++dLspJS6jv6CipkC6befWxlU+oz18mXjtLVsfWsrbiXc8H/Vn/bPtadEyeubevLU7KSUSpdtewUR+QcwBkgAfjDGTLIrluLEqf328+LIEau9P1WlT9OuPXs+3cPaqWuJ+jOKZsOb0f317pSpUSbLxRWZJi+lCjlbfjUi4g8MAJoZY2JFRAu85BOn9NvPi7SVPt9+G0aP5uTWswS2+YA/dv5BzfY1eWDVA9Rok73+/kWiyUupIsKu+wRGAbONMbEAxpjzNsVR7OR7v/28CAyE226z6v0MGgQHDxIx4GGWP7iSjzt9zNVzVxn0xSAe2/xYthMAZNzkFRmb4KQ3olTxZddhU0Ogk4jMAmKACcaYsPQmFJEngScBateuXXARFlH51m8/L1JX+mzUCEJCiGvbiU1zNrFl7hYQ6PJyFzpM7IBnKc8cL77QNnkpVQQ5bc8gIiFA1XRemuZYb3mgHXAH8LWI3GSMMWknNsYsAZYAtG7d+obX1Y1su4AaH28N7DJjBiQlwauvYp57nr3LD7Gm0UIiz0bSdFhTus/uTtlaZXO9mkLX5KVUEea0JGCM6ZHRayIyCvjWsdPfLiJJQCXggrPiUU4WGmr1+d+/H/r1g3fe4fSfHgR2/YIz289QvXV1Bn8zmFodauV5VclNXnrPgFJ5Z9evZiXQDVgvIg0BTyDcplhUXpw/DxMnwmefQZ06sGoVV1p1JWRyCL8s/YXS1Uoz4JMBNB/eHHHLn/r+eWny0l5FSl3Prm//R8BHIrIPiANGpNcUpAqxxER4/32YNs2q9Dl1KvHjJ/HT4j1sHrqQpMQkOk7tSKcpnfAsnfN2/6zkpslLexUpdSNbvvnGmDjgITvWrfJBWJg1yMvOndCtG2bhQvb/nEjI7R9z+eRlmgxuQo85PShfr3yBhZSdI/xCfSOdUjbRb77KvkuXrCP/996zKn1+9RVn63ci8IkgTm0+RdUWVRn42UDqdqlboGFl9whfexUpdSMdT0BlzRirzb9RI6sJ6NlniVy/k1VBJfmgzYf8dfgv+n3Qjyd2PFHgCQCyf99Acq+i1LRXkXJ1eiagMrdvn9XrZ+NGaNeOhP/7ka3rrrGx9ackxCbQYWIHOk3rhHdZ+46ks3uEr72KlLqRfvtV+tJU+jRLlnCgXAdWDw0h4lgEjQY0ImBuABVurpDrVeRXT53s3jdQKG6kU6qQ0W+/ul46lT7/fPB5Amds48SG5fjd5sfwkOHc1P2mPK0mP3vq5OQIXyuRKnU9KUo9M1u3bm127NhhdxjFV5pKn1dnLWDtqkh2fbgLn4o++M/0p9XIVrh55P1S0rkrMXSas+6Go/fQSf652kEnn1XoEb5S6crwJh39lagbKn0mzp3PtsTWhA7bRHx0PO3GtaPLS13wLpd/R8/53VNHj/CVyh1NAq4uMBDGjIGjRzEPDOW3Xv8geNYO/jqyhgZ9GxAwL4BKjSrl+2q1/o9ShYMmAVeVptLn+Q9WEfSfy/z+aDCVGlfiwcAHubnXzU5bvfbUsU98fDynT58mJibG7lBUPvP29qZmzZqUKJH9gym9JuBq0lT6jH5uGuvCm7Hzg914lfWi64yutH66Ne4l3J0eirbj2+PYsWP4+vpSsWLFlPGbVdFnjOHixYtERkZSr169tC/rNQHFdZU+E/v2I+z2p9nwzi/ERu6m9ejWdJ3eFZ+KPgUWjrbj2yMmJoa6detqAihmRISKFSty4ULOijFrEnAF587BpEkplT4Pv/QZwV9fIfyHMG7qeRO95vfC71Yd4TO3imJlUk0AxVNuPtfC/U1VeZOYCEuWwNSpcPUq4U9OJeh4I4688jsVGlRg6P8NpUHfBrpDyAOtTKqKOq0dVFyFhUG7djB6NNeatSXwgU949yNvTm07S8C8AEbvG03DuxtqAsgjHe845y5evEiLFi1o0aIFVatWpUaNGimP4+Li8m09ISEhiAiffvppynNhYWGICAsWLMj2co4cOUKLFi3yPE1hpYcqxU2qSp9JflXZ8chi1n8XScymI7Qc2ZJuM7tRyq+U3VEWG1qZNOcqVqzInj17AJg+fTqlS5dmwoQJ101jjMEYg5tb3o5TmzZtyrJlyxgxYgQAy5Yto3nz5nlaZnGjSaC4MAY+/xwmTICLFzk64DmCDtXlwifnqetfl94LelOlWRW7oyx2ivz9DuPGgWOHnG9atLB6oOXQkSNHGDhwIB07dmTbtm2sXLmS5s2bExERAVg78JCQED788EPOnTvHqFGjOHnyJG5ubrzzzju0a9fuhmXedNNNXLhwgfDwcCpWrMjq1au56667Ul7ftWsXo0aN4tq1azRo0ICPPvqIsmXLEhYWxuOPP06pUqW48847U6ZPSEhg0qRJbNq0iZiYGJ599llGjhyZi41UeGhzUHGwfz906QIjRnCxelO+6riYL1b6khCbxP3f3s/Dax7WBOAkyfc7eDlKaej9Dnnz66+/8vjjj7N7925q1KiR4XTPPvsskyZNYseOHXz99deZ7ojvvfdeli9fTmhoKG3btr2uD/1DDz3EvHnz2Lt3L40aNWLmzJkAPPLII7z77rts2bKFxMTElOmXLFmCn58f27dvJywsjEWLFnHy5Ml8eOf20W9qURYVZfX3X7CAmNKVCA2Yy7Z10Xh4/UWPOT1oO7YtHrozcqoiX5k0F0fszlS/fn3uuOOOLKcLCQnh0KFDKY8vXbrEtWvXKFmy5A3TDhkyhOHDh9OwYUOGDh3K2rVrAevaRExMDB07dgRgxIgRDB8+nPDwcK5du5ZyBjB8+HDWrVsHQHBwMAcOHGDZsmUAXL58mcOHD1OnTp28vXEbFZFvqrpOqkqfSafPsLvjs6w9WJ3o1VG0eLQF3Wd1p3TV0nZH6TL0fof8U6rU39er3NzcSH0za+o7nI0xbN++HU/PrMevrlGjBsYYNmzYwOLFi1OSQGY3ymbUYcIYw+LFi+nevft1zx85ciTLOAorbQ4qao4cgT594L77OO7ViCUN5vL9pvJUuqUyT+54kgH/HqAJQBULbm5ulC9fnsOHD5OUlMSKFStSXuvRoweLFi1Kebwni+saM2fOZM6cOdddaK5UqRIlS5bkp59+AuDzzz+nS5cuVKpUCW9vb7Zs2QLA0qVLU+bp1asXixcvJiHB6v116NAhrl27lvc3ayM9EygqYmJg9myYPZtLHpVZ3WwWB/bGU7a2G/f95z6aDG6i3T1VsTNnzhx69+5N7dq1adKkCbGxsQAsWrSIUaNG8fHHH5OQkIC/v/91SSGt5CaftD7//POUC8M333wzH3/8MQAff/wxI0eOpFSpUgQEBKRM/9RTT3Hy5MmU7qB+fn6sWrUqv96uLbR2UFHw44/wj38Qe/QUGxs/xdajlXHzcKfjlI60f749JUoWkZ4oqlA4cOAAjRs3tjsM5SQZfL5aO6hIOnUKxo3DfLuCPVV7sbbCY0QdiKfZ8Nvo/np3ytQoY3eESqkiTpNAYZSq0ufJ+GoEVp/OH2cNNdtVYciCXtRsW9PuCJVSxYQmgcLGUekzYv9pQqo/zv6zFfCV0gz6ogdNhzXVdn+lVL7SJFBYnDsHEycS9/kyNpfty0+e98Nf7nR+qQN3TroTz1JZd4VTSqmc0iRgt8REeP99zJSp7L16E2tKTyXysnDb0FvpMbsHZWuXtTtCpVQxpknATmFhMGoUp3eeI7DM45xJLEP1W6oz+O3e1OpQy+7olFIuQG8Ws8OlSzB6NFfa9ODb/bfwb0ZyuVR1Bnw8gJHbRmoCUMWeu7t7SvnoFi1aMHv27AynXblyJb/++mvK45deeomQkJA8xxAREcHixYtzPN/06dOZO3duus+LyHV3D8+fPx8RISdd2z/55BPGjBmT52myS88ECpIx8NlnxE+YzE8XG7HZYzxJxoOOU9vTaUonPEtru79yDSVLlszyLt9kK1eu5O6776ZJkyYAvPLKK/kSQ3ISGD16dL4sD/4uXf3iiy8CsHz58pS4CyuXSAKFYvi/ffswo0azf9NfrPYcwRVTkiYDm9DjjR6Ur1e+YGNRyiFwXCB/7vkzX5dZtUVVei/onat5J0+ezHfffYeHhwcBAQHcc889fPfdd2zYsIFXX32V//73v8ycOZO7776b++67j7p16zJs2DDWrVtHfHw8S5YsYcqUKRw5coSJEyfy9NNPExUVxYABA7h06RLx8fG8+uqrDBgwgMmTJ3P06FFatGhBz549efPNN3nzzTf5+uuviY2NZdCgQcyYMQOAWbNm8dlnn1GrVi0qV67M7bffnm78AwcOZNWqVbz44ov8/vvvlC1b9rqqpV999RWvvfYaxhj69u3LnDlzAOsO5ddff51q1arRsGFDvLy8ALhw4QJPP/10SqXSBQsWXFfaOj/YlgREpAXwHuANJACjjTHb83s9tg//56j0efatrwh068spqlO1SRUGLehN3S51nb9+pQqha9euXTcS15QpU+jZsycrVqzg4MGDiAgRERGUK1eO/v37p+z001OrVi22bNnC+PHjeeSRR9i8eTMxMTHceuutPP3003h7e7NixQrKlClDeHg47dq1o3///syePZt9+/alnJEEBwdz+PBhtm/fjjGG/v37ExoaSqlSpVi2bBm7d+8mISGBVq1aZZgEypQpQ61atdi3bx+rVq1iyJAhKaUozp49ywsvvMDOnTspX748AQEBrFy5krZt2/Lyyy+zc+dOypYti7+/Py1btgRg7NixjB8/no4dO3Ly5El69erFgQMH8vOjsPVM4A1ghjHmRxHp43jcNb9XktHwf6GT/J2bBByVPqPGTGbNn7eyh5GUquhDv9d60OLRFri56+UYZb/cHrHnVXrNQQkJCXh7ezNy5Ej69u3L3Xffna1l9e/fH7CaYqKiovD19cXX1xdvb28iIiIoVaoUU6dOJTQ0FDc3N86cOcO5c+duWE5wcDDBwcEpO+CoqCgOHz5MZGQkgwYNwsfH57r1ZeSBBx5g2bJlBAUFsWbNmpQkEBYWRteuXalcuTIADz74IKGhoQDXPT9kyBB+++03wCqZnfp6yJUrV4iMjMzWdskuO5OAAZLrHpQFzjpjJbYM/3fkCAmjn2XL6kg2uQ0hwcOT9uPa0fnFzniX1XLDSqXHw8OD7du3s2bNGpYtW8bChQtTyj5nJrnpxM3NLeX/yY8TEhJYunQpFy5cYOfOnZQoUYK6deteV5Y6mTGGKVOm8NRTT133/IIFC3J0k2a/fv2YOHEirVu3pkyZv0u75KZ0dVJSElu2bEl3nIT8Yufh6DjgTRE5BcwFpqQ3kYg8KSI7RGTHhQsXcryS5OH/UnPa8H/XrmFens6vje9lUUgj1tKDenffyjMHxhDwZoAmAKUyERUVxeXLl+nTpw8LFixIOVPw9fXN09Hv5cuX8fPzo0SJEqxbt44TJ06ku9xevXrx0UcfERUVBcCZM2c4f/48nTt3ZsWKFVy7do3IyEj+7//+L9P1lSxZkjlz5jBt2rTrnm/bti0bNmwgPDycxMREvvrqK7p06ULbtm1Zv349Fy9eJD4+nm+++SZlnoCAABYuXJjyOLsX03PCqWcCIhICVE3npWlAd2C8Mea/InI/8G+gR9oJjTFLgCVgVRHNaQzJw/9NXL6X2IQk5w3/9+OP/PnkSwSevo0T3IPfLeUZvvBubup+U/6uR6liIO01gd69ezN27FgGDBhATEwMxhjmz58PWM0rTzzxBO+88w7Lly/P8boefPBB+vXrR+vWrWnRogW33HILYA14f+edd3Lbbbdx11138eabb3LgwAHat28PQOnSpfniiy9o1aoVQ4YMoUWLFtSpU4dOnTpluc4HHnjghueqVavG66+/jr+/P8YY+vTpw4ABAwCre2n79u2pVq0arVq1ShnS8p133uGZZ56hWbNmJCQk0LlzZ957770cb4PM2FZKWkQuA+WMMUasc6HLxphMy2LmtpR0cu8gpwz/d+oUV59+nrX/i2EXrShZpgTd5vSi1chWuHlou78qfLSUdPFWlEpJnwW6AOuBbsBhZ63IKcP/xcWROHc+26b/SGh8B+LdvGj7TGu6zOhGyfLOa79TSqn8ZGcSeAJ4W0Q8gBjgSRtjyRGzfj2/PTyL4FON+Qt/GnSrScDiAVRqVMnu0JRSKkdsSwLGmE1A+p1tC6tz5zj/xFSC/i+e3+lIpZpePPjBfdzc+2a7I1NKqVxxiTuG8ywxkeh577Lun2vYGdcML2+h98xutB7bAfcS7nZHp5RSuaZJIAuJP20lbMg8NpyuT6w0p/WDDen69kB8KvrYHZpSSuWZJoGMXLrE4YemE/y/RMK5jZua+tDry4fxu62K3ZEppVS+0SSQljGEz/6Q4Ok/cTiuLhXKJTL0/QE0GNxch3ZUKp+4u7vTtGlT4uPj8fDwYMSIEYwbNw43t4y7VR8/fpyffvqJYcOGFWCkxZ8mgVSubd7FhvsXEXa2BiXca9Bz/K20nT0Id09t91euyxlVeFPXDjp//jzDhg3j8uXLKVU703P8+HG+/PJLTQL5TJMAkHTpMjsHz2bdmiRiqEnLLr50W/YUpar62h2aUrYqiCq8fn5+LFmyhDvuuIPp06dz4sQJhg8fztWrVwFYuHAhHTp0YPLkyRw4cIAWLVowYsQIBg0alO50KmdcOwkYw9GXPiVo9m4uJFSgbrVYen01jKpdGtkdmVKFQkFV4b3ppptISkri/Pnz+Pn5sXr1ary9vTl8+DBDhw5lx44dzJ49m7lz5/L9998DEB0dne50KmdcNglcDNlF8LBP+e1CBcp7enH/6y245YX+2u6vVCoFWYU3uYRNfHw8Y8aMYc+ePbi7u6eUVU4ru9OpzLlcEoj5M4LQexawbUsSHpSme7+StPvqBTxKaYVPpdJKrsKbOhE4owrv77//jru7O35+fsyYMYMqVarw888/k5SUhLd3+r/N+fPnZ2s6lTmXSQJJiUnsHv85axcdIDrJmxY3RdL92zGUbl7f7tCUKrQKogpv8hCKY8aMQUS4fPkyNWvWxM3NjU8//TSlomba0s8ZTadyxiWSgElK4tOqkzkZXoraJaPo9VZ7qj89wO6wlCr0fLw86NG4CqGT/PO1Cm9yKenkLqLDhw/nueeeA2D06NHce++9fPPNN/j7+1OqVCkAmjVrhoeHB82bN+eRRx7JcDqVM7aVks6N3JaSBtjR5yVKlvOiyccTkFSjDynlarSUdPFWlEpJF6jW/3vF7hCUUqrQ0VFPlFLKhWkSUMoFFaVmYJV9uflcNQko5WK8vb25ePGiJoJixhjDxYsXc9xV1mWuCSilLDVr1uT06dNcuHDB7lBUPvP29qZmzZo5mkeTgFIupkSJEtSrV8/uMFQhoc1BSinlwjQJKKWUC9MkoJRSLqxI3TEsIheAE3bHkQOVgHC7g7CZbgOLbgfdBsns2A7hxpje6b1QpJJAUSMiO4wxre2Ow066DSy6HXQbJCts20Gbg5RSyoVpElBKKRemScC5ltgdQCGg28Ci20G3QbJCtR30moBSSrkwPRNQSikXpklAKaVcmCYBJxORN0XkoIjsFZEVIlLO7pgKmogMFpH9IpIkIoWma1xBEJHeInJIRI6IyGS747GDiHwkIudFZJ/dsdhFRGqJyDoROeD4LYy1O6ZkmgScbzVwmzGmGfAbMMXmeOywD7gHCLU7kIIkIu7AIuAuoAkwVESa2BuVLT4B0r1RyYUkAM8bYxoD7YBnCst3QZOAkxljgo0xCY6HW4Gc1XktBowxB4wxh+yOwwZtgCPGmN+NMXHAMmCAzTEVOGNMKPCX3XHYyRjzhzFml+P/kcABoIa9UVk0CRSsx4Af7Q5CFZgawKlUj09TSH74yj4iUhdoCWyzNxKLjieQD0QkBKiazkvTjDGrHNNMwzolXFqQsRWU7GwDFyTpPKd9sl2YiJQG/guMM8ZcsTse0CSQL4wxPTJ7XURGAHcD3U0xvTEjq23gok4DtVI9rgmctSkWZTMRKYGVAJYaY761O55k2hzkZCLSG3gB6G+MibY7HlWgwoAGIlJPRDyBB4DvbI5J2UBEBPg3cMAY85bd8aSmScD5FgK+wGoR2SMi79kdUEETkUEichpoD/wgIkF2x1QQHB0CxgBBWBcCvzbG7Lc3qoInIl8BW4BGInJaRB63OyYb3AkMB7o59gN7RKSP3UGBlo1QSimXpmcCSinlwjQJKKWUC9MkoJRSLkyTgFJKuTBNAkop5cI0CSillAvTJKCUUi5Mk4BSSrkwTQJK5ZKIVBaRP0TkpVTPNRORGBG5z87YlMouvWNYqTwQiAolzAAAAOJJREFUkV7A/wFdgD3ADmC7MeZRWwNTKps0CSiVRyKyAOgPbAA6AS2MMVH2RqVU9mgSUCqPRMQL+BloAHQwxhSKwUKUyg69JqBU3tXFGjfAADfZG4pSOaNnAkrlgWOgkC3AYazhAqcDzYwxJ+2MS6ns0iSgVB6IyGxgGNAMuIw1hnRJwN8Yk2RnbEplhzYHKZVLItIFeB542BgT4Rg69BGgMdZockoVenomoJRSLkzPBJRSyoVpElBKKRemSUAppVyYJgGllHJhmgSUUsqFaRJQSikXpklAKaVcmCYBpZRyYf8PtOVT54GYFWgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate data\n",
    "np.random.seed(123)\n",
    "N = 20\n",
    "beta0 = -4\n",
    "beta1 = 2\n",
    "x = np.random.randn(N)\n",
    "e = np.random.randn(N)\n",
    "y = beta0 + beta1*x + e\n",
    "true_x = np.linspace(min(x), max(x), 100)\n",
    "true_y = beta0 + beta1*true_x\n",
    "\n",
    "# estimate model \n",
    "beta1_hat = sum((x - np.mean(x))*(y - np.mean(y)))/sum((x - np.mean(x))**2)\n",
    "beta0_hat = np.mean(y) - beta1_hat*np.mean(x)\n",
    "fit_y = beta0_hat + beta1_hat*true_x\n",
    "\n",
    "# plot\n",
    "fig, ax = plt.subplots()\n",
    "sns.scatterplot(x, y, s = 40, label = 'Data')\n",
    "sns.lineplot(true_x, true_y, color = 'red', label = 'True Model')\n",
    "sns.lineplot(true_x, fit_y, color = 'purple', label = 'Estimated Model')\n",
    "ax.set_xlabel('x', fontsize = 14)\n",
    "ax.set_title(f\"Linear Regression for y = {beta0} + {beta1}x\", fontsize = 16)\n",
    "ax.set_ylabel('y', fontsize=14, rotation=0, labelpad=10)\n",
    "ax.legend(loc = 4)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<u>Extensions of Ordinary Linear Regression</u>**\n",
    "\n",
    "- Regularized\n",
    "- Bayesian\n",
    "- GLMs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
